function [spikes, npmin, vdmin, vdmean] = srf_FindSpikes(hfile, cutoff)
% SRF::FindSpikes  - find spiked vertices (1-based indices)
%
% FORMAT:       [spikes, npmin, vdmean] = srf.FindSpikes([cutoff]);
%
% Input fields:
%
%       cutoff      cutoff value for critical sqrt(abs(n*n')),
%                   default: 1/3
%
% Output fields:
%
%       spikes      Vx1 vertex list of found spikes
%       npmin       Vx1 minimal product list per vertex
%       vdmin       Vx1 min vertex-to-neighbors distance
%       vdmean      Vx1 mean vertex-to-neighbors distance
%
% Note: the srf object remains unchanged, often it is more useful
%       to apply e.g. coloring to an unsmoothed mesh

% Version:  v0.7b
% Build:    7090215
% Date:     Sep-02 2007, 3:50 PM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% check arguments
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'srf')
    error( ...
        'BVQXfile:BadArguments', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if nargin < 2 || ...
   ~isa(cutoff, 'double') || ...
    numel(cutoff) ~= 1 || ...
    isinf(cutoff) || ...
    isnan(cutoff) || ...
    cutoff <= 0 || ...
    cutoff >= 1
    cutoff = 1/3;
else
    cutoff = real(cutoff);
end

% get vertices, normals, and neighbors
vx = bc.VertexCoordinate;
no = bc.VertexNormal;
nb = bc.Neighbors;
nv = size(no, 1);

% initialize array
npmin = ones(nv, 1);
vdmin = zeros(nv, 1);
vdmean = zeros(nv, 1);

% iterate over vertex normals
for nc = 1:nv
    
    % get minimal product of this vertex' neighbors with normal
    tnb = nb{nc, 2};
    npmin(nc) = min(sqrt(abs(no(nc, :) * no(tnb, :)')));
    
    % calc mean vertex-neighbor distance
    vd = sqrt(sum( ...
        (repmat(vx(nc, :), [numel(tnb), 1]) - vx(tnb, :)) .^ 2, 2));
    vdmin(nc) = min(vd);
    vdmean(nc) = mean(vd);
    
end

% find spikes
spikes = find( ...
    (npmin .* (vdmean / mean(vdmean)) .* (vdmin / mean(vdmin))) <= cutoff);
